/*!  \file  GeSVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/28/2013 12:19:16 PM
 *   \copyright GNU Public License.
 */
#include "FullSVD.h"
#include "Matrix.h"
#include "GeSVD.h"
#include "setThreads.h"
extern "C" {
void dgesvd_  ( //computes SVD MxN matrix A, A = U * SIGMA * transpose(V)
const char   * jobu , // 'A': All columns, 'S'  the first min(m,n) of U/V are returned in array U/V
const char   * jobvt, // 'O':  the first min(m,n) rows of U/V^T, 'N' : no singular vectors
const int    *     m, // Rows of the input matrix A
const int    *     n, // Columns of the input matrix A
      double *     A, // (input/output) On entry, the M-by-N matrix A., at output content is destroyed
const int    *   ldA, // leading dimension of the array A, should be M
      double *     s, // double[min(M,N)] containg singular values
      double *     u, // left singular vectors, stored columnwise
const int    *   ldU, // 
      double *    vt, // right singular vectors, stored rowwise
const int    *  ldvt, // 
      double *  work, // working array of size lwork
const int    * lwork, // The dimension of the array WORK. LWORK = -1,then a workspace query is assumed
int    *  info  //  0  :  successful exit.
                // < 0 :  if INFO = -i, the i-th argument had an illegal value.
                // > 0:  if DBDSQR did not converge, INFO specifies how many
                // superdiagonals of an intermediate bidiagonal form B
                // did not converge to zero. See the description of WORK
                // above for details. 
);
}
GeSVD::GeSVD(): FullSVD( "Lapack dgesvd" ) {
	aColMajor = true; aRowMajor = false;
	setCompileTime(__DATE__, __TIME__);	
	canBeParallel = true;
}
bool GeSVD::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
void GeSVD::compute(Matrix & A, Matrix & U, Matrix & Sigma, Matrix & VT ) {
	reset();
	int m = A.rows(), n = A.cols(), lda = A.lda(), ldu = U.lda(), ldvt = VT.lda();
	int info = 0, lwork = -1;
	double* work = 0, wkopt = 0;
	char job = 'A';
	dgesvd_(&job, &job, &m, &n, A.getModData(), &lda,	Sigma.getModData(), 
			U.getModData(), &ldu, VT.getModData(), &ldvt, &wkopt, &lwork, &info );
	lwork = (int) wkopt;
	allocate(work, lwork);
	if (work) { /*  Compute SVD */
		dgesvd_( &job, &job, &m, &n, A.getModData(), &lda,	Sigma.getModData(), 
					U.getModData(), &ldu, VT.getModData(), &ldvt, work, &lwork,
					 &info );
		delocate(work);	
	};
};

REGISTER(GeSVD);
